#include <complex.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <fftw3.h>

#define pi 2*acos(0)

int main(int argc, char** argv)
{
  fftw_complex wrk1, wrk2;
  FILE *fh;
  char param[80], value[80], line[80];
  int j, n, ic_switch, kp, kz, w_flag;
  double l, ap, az, x, dx, az2, kz2;

  w_flag = 0;
  if (argc != 2) 
  {
    printf("Usage:\t%s input_file\n",argv[0]);
    exit(1);
  }
  fh = fopen(argv[1],"r");
  if (fh == NULL) 
  {
    printf("File unreadable\n");
    exit(1);	
  }
  while (fgets(line, 80, fh) != NULL)
  {
    sscanf(line, "%s\t%s", param, value);
      if (strcmp(param,"#gridpoints") == 0) n = atoi(value);
      if (strcmp(param,"#box_length") == 0) l = atof(value);
      if (strcmp(param,"#height_psi") == 0) ap = atof(value);
      if (strcmp(param,"#height_z") == 0) az = atof(value);
      if (strcmp(param,"#height_z2") == 0) az2 = atof(value);
      if (strcmp(param,"#wavenum_psi") == 0) kp = atoi(value);
      if (strcmp(param,"#wavenum_z") == 0) kz = atoi(value);
      if (strcmp(param,"#wavenum_z2") == 0) kz2 = atoi(value);
      if (strcmp(param,"#ic_type") == 0) ic_switch = atoi(value);
      if (strcmp(param,"#1.") == 0) w_flag = 1;
  }
  fclose(fh);    
  printf("Read input file\nGridpoints\t%d\nBox length\t%f\n",n,l);
  printf("IC type   \t%d\nHeight of psi\t%f\nHeight of z\t%f\n",ic_switch,ap,az);
  printf("Wavenum psi\t%d\nWavenum z\t%d\n", kp, kz);
  if (w_flag == 0) {
    printf("Generating IC\n");
    dx = (l/n);    
    fh = fopen(argv[1],"a"); 
    fprintf(fh,"## IC generated by %s\n#1. x\t\t\t2. Psi_re\t\t3.Psi_im",argv[0]);
    fprintf(fh,"\t\t4. z_re\t\t\t5. z_im\n\n");
    
    for (j = 0; j < n; j++ )
    {
      x = -l/2 + (j+0.5)*dx;
      wrk1 = ap*cexp(1.I*(2*pi/l)*kp*x);
      wrk2 = az*cexp(1.I*(2*pi/l)*kz*x) + az2*cexp(1.I*(2*pi/l)*kz2*x);
      fprintf(fh,"%.15e\t",x);
      fprintf(fh,"%.15e\t%.15e\t",creal(wrk1),cimag(wrk1));
      fprintf(fh,"%.15e\t%.15e\n",creal(wrk2),cimag(wrk2)); 
    }
    fclose(fh);
  } else {
    printf("IC generated.\nNothing to do.\n");
  }
  	
  return 0;
}

